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Abstract 

A computer code for quasiparticle random phase approximation-QRPA and projected quasipar- 
ticle random phase approximation-PQRPA models of nuclear structure is explained in details. The 
residual interaction is approximated by a simple 5-force. An important application of the code con- 
sists in evaluating nuclear matrix elements involved in neutrino-nucleus reactions. As an example, 
cross section for 56 Fe and 12 C are calculated and the code output is explained. The application to 
other nuclei and the description of other nuclear and weak decay processes is also discussed. 



Program summary 

Title of program: QRAP (Quasiparticle RAndom Phase approximation) 

Computers: The code has been created on an PC, but also runs on UNIX or LINUX machines. 
Operating systems: WINDOWS or UNIX 
Program language used: Fortran-77 

Memory required to execute with typical data: 16 Mbytes of RAM memory and 2 MB of hard disk 
space 

No. of lines in distributed program, including test data, etc.: ~ 8,000 
No. of bytes in distributed program, including test data, etc.: ~ 256 kB 
Distribution format: tar.gz 

Keywords: QRPA; Projected QRPA; semileptonic processes. 

Nature of physical problem: The program calculates neutrino- and antineutrino-nucleus cross 
sections as a function of the incident neutrino energy, and muon capture rates, using the QRPA or 
PQRPA as nuclear structure models. 

Method of solution: The QRPA, or PQRPA, equations are solved in a self-consistent way for 
even-even nuclei. The nuclear matrix elements for the neutrino- nucleus interaction are treated as 
the beta inverse reaction of odd-odd nuclei as function of the transfer momentum. 

Typical running time: w 5 min on a 3 GHz processor for Data set 1. 

PACS numbers: 21.60.Jz, 25.30.Pt, 26.30.Jk 



Long Write-Up 

I. INTRODUCTION 

The new age of the physics beyond the standard model 
of electroweak interaction has as one of the most promis- 
ing pathways the search of neutrino oscillations. Sev- 
eral experimental efforts are oriented to find the neutrino 
masses and the related oscillations involvingatmospheric, 
solar, reactor and accelerator neutrinos [lHy]- Since neu- 
trinos interact so weakly with matter, they bring infor- 
mation on the dynamics of supernova collapse and pos- 
terior explosion as well as on the synthesis of heavy nu- 
clei [10|. 
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The detection signal of neutrinos is measured trough 
the weak interaction of incoming neutrinos with the nu- 
clei present in, e.g., a liquid scintillator detector, as well 
as with the surrounding blockhouse detector-shield. The 
flux-averaged i/-nucleus cross sections are the measured 
observables. Recently, Ref. [|| has studied the effect of 
neutrino oscillations on the expected supernova neutrino 
signal with the LVD detector, through their interactions 
with protons and carbon nuclei in a liquid scintillator and 
with iron nuclei in the support structure. 

Charged and neutral i/ e -nucleus cross sections on 12 C 
(liquid scintillator) as well as on 56 Fe (detector surround- 
ing shield) were measured by the KARMEN Collabora- 
tion H,[lJ|- Other experiments such as LAMPF [O, OJ] 
and LSND 03] have also used 12 C to search for neu- 
trino oscillations and to measure neutrino-nucleus cross 
sections. Furthermore, future experiments will use 12 C as 
liquid scintillator, such as in the spallation neutron source 
i{ibblM!ri^i&fc.&Iakctti^ RNL) , 
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or in the LVD (Large Volume Detector) experiment Q . 

On the other hand, the cross sections v e {v e )— 56 Fe are 
important to test the ability of nuclear models in explain- 
ing reactions on nuclei with masses around iron, which 
play an important role in supernova collapse [16]. The 
iron is used as material detector in experiments on neu- 
trino oscillations such as MINOS [171. whereas future ex- 
periments, such as SNS at ORNL [15| plan to use the 
same material. 

There have been great efforts on nuclear structure 
models to describe consistently semileptonic weak pro- 
cesses with 12 C such as RPA-like models. A brief sum- 
mary on the different models employed for 12 C is sketched 
in Ref. [l|- 

The puzzle with the Random Phase Approximation 
(RPA) and the quasiparticle RPA (QRPA), when applied 
to the weak observables in the triad { 12 B, 12 C, 12 N}, is 
well known. That is, to get agreement with data for 
the ground state triplet T = 1 {P ± -decays, /i-capture, 
and the exclusive 12 C(^ e , e~) 12 N reaction) the contin- 
uum RPA (CRPA) calculations of Kolbe, Langanke, and 
Krewald [19( needed to be rescaled by a reduction fac- 
tor = 4. The reason for such a large discrepancy is very 
simple: within the RPA the transitions 12 C— ;> 12 N(1+) 
and 12 C— s- 12 B(l+) are engendered mostly by the particle- 
hole excitation p 3 / 2 — > Pi/2> what is physically incorrect. 
In fact, since late 1980's we know from several hadronic 
charge-exchange reaction measurements, and the consec- 
utive Shell Model (SM) calculations, that the excitations 
P3/2 ~> P3/2, Pi/2 -> Pi/2, and p 1/2 -> p 3/2 partici- 
pate quite significantly in these processes (see, for in- 
stance, [2(| Table I]). It is the involvement of these con- 
figurations that brings about the necessary quenching of 
the Gamow- Teller (GT) resonances and /3-decay rates. 
To make them come into play it is mandatory to open the 
P3/2 shell by means of pairing correlations, which is done 
within both the SM and the QRPA. But, a new problem 
emerges in the application of the QRPA to 12 C, as first 
observed by Volpe et al. 21] who noted that within this 
approach the lowest state in 12 N irremediable turned out 
not to be the most collective one. As a consequence the 
QRPA also fails in accounting for the exclusive processes 
to the isospin triplet T = 1. Soon after it was shown 
[13-113] that the origin of this difficulty arises from the 
degeneracy among the Px/2 an d Pz/2 quasiparticle ener- 
gies (both for protons and neutrons), which is inherent 
to the non-conservation of particle number. Therefore, 
for a physically sound description of the weak processes 
among the A = 12 iso-triplet it is imperative to use the 
SM or the number projected QRPA (PQRPA). 

The QRAP code is based on Refs. [H-|13, where a 
new formalism for neutrino-nucleus scattering has been 
developed, and the PQRPA is used as the nuclear model 
framework. The residual interaction was done with the 
simple 5-force, which has been used extensively in the lit- 
erature to describe the single and double beta decays [25l — 

Before proceeding we address briefly on the genesis of 



the QRPA and PQRPA in a manner appropriate in the 
present context. Although this is not a topic of central 
interest for the application-oriented computer code, it 
belongs to the physics background. The neutron-proton 
QRPA was developed in 1967 by Hableib and Soren- 
son [3l| in order to account for the hindrance of the 
allowed /3-transitions. Almost 20 years later, when Vo- 
gel and Zirnbauer [32| and Cha [33J discovered the im- 
portance of the particle-particle force in the S = 1, T 
= channel, the QRPA became to be the most fre- 
quently used nuclear structure method for evaluating 
double beta (/?/3) rates. It was quickly realized, how- 
ever, that a small change in the particle-particle interac- 
tion strength caused a large change in the lifetimes and 
eventually the breakdown (called a collapse) of the entire 
method. Later on several modifications of the QRPA 
were proposed to make it more reliable. One of these 
was the charge-exchange PQRPA, which has been for- 
mulated to evade the disadvantages inherent in the non- 
conservation of particle number, and was derived from 
the time-dependent variational principle 29]. But, the 
PQRPA did not yield substantially different result from 
the plain QRPA, and was unable to avoid the collapse in 
the study the two-neutrino /3/3-decay in 76 Ge. As a mat- 
ter of fact, the problem of the QRPA collapse has not yet 
been settled down, in spite of enormous effort invested for 
this purpose by many nuclear physicists (compare, for in- 
stance, Fig. 1 from Ref. [2!| with Fig. 5 from a recent 
work of Yousef et al. [34|). 

However, the PQRPA turned out to be quite important 
for the description of relatively light nuclei such as 12 C. 
For example, the employment of PQRPA for the inclusive 
12 C(i/ e , e~) 12 N cross section, instead of the continuum 
RPA (CRPA) used by the LSND collaboration in the 
analysis of — > v e oscillations of the 1993-1995 data 
sample, leads to an increased oscillation probability [!ij . 

The PQRPA was recently also used to calculate the 
56 Fe(i/ e , e~) 56 Co cross section (35|. A comparison be- 
tween the QRPA and PQRPA for the same interaction 
and employing the same model space shows that the pro- 
jection procedure could be important for medium mass 
nuclei. Moreover, several approximations such as: i) Hy- 
brid Model (HM) [Hj], ii) QRPA with Skyrme interac- 
tion [371, hi) relativistic QRPA (RQRPA) [Hj], and iv) 
QRPA and PQRPA with the <5-force [H| yield different 
results for the neutrino cross section as a function of the 
neutrino energy. It is a hard task to find the origin for 
the differences, mainly because these models are not us- 
ing the same interaction and/or the same single-particle 
configuration space, carrying different types of correla- 
tions in each case. 

The cross sections for charged- and neutral-current 
neutrino-induced reactions on the iron isotopes 52_60 Fe 
were also evaluated within the HM for various supernova 
neutrino spectra [3{|- Here, large-scale SM calculations 
were used for the GT-like contributions, while transitions 
for other multipoles are based on the RPA. More pre- 
cisely, the authors scale the SM cross sections using the 
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ratios obtained from the RPA calculations with and with- 
out this dependence of the multipole operator. The rea- 
son for such a procedure is twofold: i) the limitation of 
the SM to account for momentum-transfer dependence of 
the GT operator, and ii) the lack of pairing correlations 
in the RPA. It should be also mentioned that SM calcu- 
lations of inelastic neutral-current neutrino-nucleus cross 
sections in medium-mass nuclei, present in supernova en- 
vironment, have been constrained by the highly precise 
data on the magnetic dipole strength distributions for 
the nuclei 50 Ti, 52 Cr, and 54 Fe, which are dominated 
by spin-isospin flipping (GT-like) contributions [40| . In 
spite of the agreement between data and calculations it 
was necessary to consider also here the effects of finite 
momentum transfer what was done via the RPA. Briefly, 
the HM is neither fish nor fowl, and a comparison of the 
results from Refs. [39l . Eoj with self-consistent calcula- 
tions, such as the QRPA, PQRPA and RQRPA, could 
be enlightening. 

This brief introduction shows: 1) the importance of 
neutrino-nucleus cross sections for astrophysical purposes 
and, 2) that these cross sections are strongly correlated 
with the nuclear structure model employed. The QRAP 
code, with a simple residual interaction, is able to access 
the sources of these problems and it can calculate several 
weak interaction processes mentioned above. Needless 
to stress that this code can be easily adapted for the 
evaluation of /3/3-decays. 

The write-up is organized as follows. In section II we 
make a short survey of the theoretical description of weak 
interaction processes, with emphasis on the formulation 
implemented in this numerical code. In sections III and 
IV we describe the QRPA, and PQRPA formalisms, mak- 
ing explicit the differences among them. In section V we 
show how the code is organized, how to make an input 
and how to understand the output. Section VI explains 
the role of each subroutine of the code. Finally, section 
VII proposes a few cases to practice with the code. 



II. WEAK INTERACTING PROCESSES 

In this section we give a brief summary of the main 
formulae developed in Ref. [HI, [23[ for: 

• neutrino scattering (NS) 

v t + (Z,A) -> (Z + l,N-l)+r, 

• antineutrino scattering (AS) 

P e + (Z,A) -> (Z -1,N + l)+£+, 

• muon capture (MC) rate 

fj,-+(Z,A) -> (Z-l,jV-l) + i/ M , 

where I = e, /z. The comparison with other for- 
malisms [4l| - |43| can be found is in just mention works. 



The weak Hamiltonian is expressed in the form 

G k 



(2.1) 



where G = (3.04545 ± 0.00006) xlfr 12 is the Fermi cou- 
pling constant (in natural units), 



J a = (J,iJ$) 

= ili 



(2.2) 



2JV1 me 



is the hadronic current operator 1 , and 

l a ((l,E v ) = = -iu se (p, £ , < .)7 Q (1 + 7 5 )m Si ,, 

(2.3) 

is the plane wave approximation for the matrix element 
of the leptonic current in the case of neutrino reactions, 
with pe = {p, iEg} and q v = {q, iE v } being, respectively, 
the lepton and the neutrino momenta. 

For the sake of convenience we will use spherical co- 
ordinates (to = —1,0,-1-1) for the three- vectors, and the 
Walecka's notation [42|], with the Euclidean metric, for 
four-vectors, i.e., x — {x, £4 = ixq,}. The only differ- 
ence is that we substitute Walecka's indices (0, 3) by our 
indices (0,0), i.e. we use the index for the temporal 
component and the index for the third spherical com- 
ponent. 

The quantity 



k = Pi-P f = {k,ik^}, 



(2.4) 



is the momentum transfer, where Pj and Pt are momenta 
of the initial and final nucleus, M is the nucleon mass, 
THe is the mass of the charged lepton, and g v , g A , g M 
and <?p are, respectively, the vector, axial-vector, weak- 
magnetism and pseudoscalar effective dimensionless cou- 
pling constants. Their numerical values are: 



9v 



1; g A = 1.26; 



9m — K p K n — 3.70; g P — 9^- 



2Mm/ 



(2.5) 



In the numerical calculations we use an effective axial- 
vector coupling g A = 1 [iij . 

The finite nuclear size (FNS) effect is incorporated via 
the dipole form factor with a cutoff A = 850 MeV, i.e., 



9^9 



A 2 



A 2 + fc 2 



(2.6) 



To use (|2.1j) with the non-relativistic nuclear wave 
functions, the Foldy-Wouthuysen transformation has to 



To avoid confusion, we will be using roman fonts (M,m) for 
masses and math italic fonts (M,m) for azimuthal quantum num- 
bers. 



4 



be performed on the hadronic current (I2.2|) . When the 
velocity dependent terms are included this yields (45j : 

J$ = 9v + (g A + 9pi)<r ■ k - g A cr ■ v, 
J = -g A a - ig w a xk-g v k + g P2 ((T ■k)k + g v v, 

(2.7) 

where k = k/ft, k = |k|, and v = — iV/M is the veloc- 
ity operator, acting on the nuclear wave functions. The 
following short notation 

,9v = 5 V 2M' ^ A = ^ A 2M' ^ w = v 2M' 

5pl = 5p 2MnV 5 ^ 9p 2M^' (2 - 8) 

has also been introduced. 

In performing the multipole expansion of the nuclear 
operators 

O Q = (O,O ) = J a e~ ikr , (2.9) 

it is convenient: 

1) to take the momentum k to be along the z axis, i.e., 



-ik-r 



(2.12) 



£rV442L + l)j L (p)*lo(*), 

L 

= £rV47r(2J + l)jj(p)y, (f), (2.10) 
J 

where p — nr, and 

2) to introduce the operators Q j, defined as 

Q = (O, Of)) = i _J V2JTT0aj. (2.11) 

J 

Thus, 

Oflj = jj(p)Yjo(f)J , 
roJ = X) iJ " LF ™uiL(/»)[il(*)®J]j, 

L 

where the geometrical factors 

are listed in Table I of Ref . [23( . 
Explicitly, from (|2~71) 

O 0J - 5v^WT + ig A Mj + i(g A + g P1 )Mtj (2.14) 

O m j = j(<5mOff P2 - .9a + mg w )M^ 

+ 5v ^^j - S m0 g v M]. (2.15) 

The elementary operators are given by 

M] = jj(p)lj(r), 

= ^i J - L - 1 -FmuiL(p)[il(f)®o-] J , (2.16) 

L>0 

L>0 



Here we make use of the conserved vector current 
(CVC). From (|233)> . (|2~T51 . and [H, Eq. (10.45) and 
(9.7)] 



k o v = ko v = k $ o; 



which yields 



QvMli - g v M] = -^g v M) 



(2.17) 



(2.18) 



Therefore, from (|2.15l) 



O mJ = i(S m0 g P2 - g A + mg w )M^ 

+ 2\m\g v M l mJ +S mQ '^g v M]. (2.19) 



The elementary operators Mj , Mj, M-qj and VWqj are 
real, but M±u and -M^u are not, and it is convenient to 
put in evidence their real and imaginary parts, expressing 
them as 



M±u = M?j±iM[ 



(2.20) 



with AJfj, and M[j arising, respectively, from the terms 
in (gUm with L = J ± 1, and L = J. Note that F± w = 
T1/V2. 

It is also convenient to separate the elementary opera- 
tors into: 

• natural parity (NP), (it = (-) J ): M] , Mfj 1 , and 
M^f, and 

• unnatural parity (UP), (tt = (-) J+1 ): Mj, M^j, 

The operators Q j = (O0j,O m j) can be express as a 
sum of real and imaginary operators, i.e., Q j = O^j + 
iO' aJ , with O^j (O^j) being a NP (UP) operator. This is 
a very important finding because it implies that O^j and 
Oqj do not contribute simultaneously, and, therefore, one 
always can deal only with real operators. 

In summary, natural and unnatural parity operators 
are, respectively: 



Ci" 



0J 



9vM] 



O *j = -^g v M], 



o r 



m^OJ 



(mg A -gJMfj I + g v M^f, (2.21) 



and 



O0j = g A Mj + (g A + g P1 )Mh, 

Ooj = (5p2-5a)Xoj, 

{-g A + mg^)Mif + g v M v ^. (2.22) 
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A. Neutrino-nucleus cross section 

For the neutrino-nucleus reaction, the momentum 
transfer is k = pi — g„, and the corresponding cross sec- 
tion reads 



<E t ,J f ) 



\Pe\E, 
2vr 



■F(Z±l,E e )J d{cos9)%(q,Jf), 

(2.23) 



where F(Z ± l,Ei) is the Fermi function (Z + 1, for 
neutrino, and Z — 1, for antineutrino) , 9 = q • p is the 
angle between the incident neutrino and ejected lepton, 
and the transition amplitude is 



Tct(k ' J/) = 2j^TT^ E KWI^IW 



(2.24) 



After expressing the spatial part of the lepton traces C a p 
in spherical coordinates, and applying the Wigner-Eckart 
theorem, one can cast the transition amplitude in the 
compact form [23[ 



%(K,Jf) 



AttG' 
2J,.+ 



yED^/IIOwII^I 2 ^ 



+ Yl K^/IIOmjll^)! 2 /:™ (2.25) 

m=0±l 

- 2»(| < J> 1 1 O flJ 1 1 Ji>< J> 1 1 Ooj 1 1 Ji» J C 0O ] • 

The explicit expressions for the traces £<$ = £00, £ m = 
£ mm , and £00 are [23] 



£00 — 1 
£00 = 
Co = 1 



|p| cos 9 

qo_ po\ 

E v Ei ) ' 

, 2q p |p|cos6> 



EfE,, 



E, 



with 



% = fc • q 



Po = ■ p 



EJMcoae-Eu) 



|p|(|p| - £„cos< 



(2.27) 



being the z-components of the neutrino and lepton mo- 
menta, and Si = ±1 for NS and AS, respectively. 



B. /^-capture rates 



reaction amplitude, by keeping in mind that: i) the roles 
of p and q are interchanged within the matrix elements 
of the leptonic current, which makes that in (|2.26j) Si —> 
— 1, ii) the momentum transfer turns out to be k = q — p, 
and therefore the signs on the right-hand sides of (qo,Po) 
have to be changed, and iii) the threshold values (p — ► 
: q —> k, &0 — s- E v — mi) must be used for the lepton 
traces. All this yields qo = E v , po = 0, and 



£n 



C 



00 



£ = 1, £1 = 0, £_i 



(2.28) 



Instead of summing over the initial lepton spins sj>, as 
done in (|2.24[) . one has now to average over the same 
quantum number. We get 



A(J/) 



E 
2tt 



^is\ 2 T MC (J f ), 



(2.29) 



where (f>ig is the muonic bound state wave function eval- 
uated at the origin, and E v = m M — (M„ — M p ) — Eg — 
Ef + Ei, where Eg is the binding energy of the muon in 
the IS" orbit. Thus from (1235)) and (|2~2%)) 



%(K,Jf) 



AnG 2 



2Ji + 1 EDWIo 0J 

2|(J / ||0_ 1J ||J,)| 2 ]. 



Ooj||^i, 



(2.30) 



In the case of MC it is convenient to rewrite the effec- 
tive coupling constants (|2.8I) as 



fh 
9* 



El - 

9v 2M' 9a 



9t 



El- 

2M' 



(,9v + 5m)^; g P =g P ^, (2.31) 



where g p = g P2 - g pl . 2 

Thus, natural and unnatural parity operators are now, 
respectively: 



Cl R — (~) R 
U 0J U 0J 

U -1J 



9v)M) 



E v 



j > 



-{gA + g^Mff + gyM^f, (2.32) 



and 



O^j - O 0J = g A M] + {g A +g A -g P )M^ 

Oiu = -(g A + g w )Mlf~g v M^. (2.33) 



III. NUCLEAR STRUCTURE CALCULATION 
A. PQRPA 

The PQRPA for charge-exchange excitations was de- 
rived from the time-dependent variational principle in 



The muon capture transition amplitude T MC (Jf) can 
be derived from the result (|2.25[) for the neutrino-nucleus 



2 Note that there is a misprint in Eq. (2.41) of Ref. [gjj. Also in 
Eq. (2.42) of the same reference g pl should read g p . 
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Ref. [29| . In the same reference is also described in details 
the projected Barden-Cooper-Schiffer (PBCS) approxi- 
mation. Basically one employs the number projection 
operators Pj\f on the \BCS) state. That is: Pq = PzPn 
for a ground state with Z protons and N neutrons, and 
P^ = Pz+hPn—h, with fi = ±1, for excited states in 
nuclei with Z + protons and N — fj, neutrons. In this 
section we give a brief description of both the PBCS and 
PQRPA approximations. 

The PBCS gap equations are 



2e k u k v k - A k (u 



= 0, 



(3.1) 



where 



Uk'Vk'yjyKkk k ; 0) 



2^ (2j fc + l)V2 



I z 



(3.2) 



are the pairing gaps, and 

I z -\k) , ^ (2 Jfc , + l) 1/2 



e k = efe- 



E 



I z ' ^ (2 Jfc + l)V2^' 



xF(kkk'k';0) 



I z 



Ae fc , (3.3) 



are the dressed single-particle energies, where 



1 f dz 



I (hk 2 --k n ) = —j-^ l a kl ---a kn 



r -l — „,2 i „2„,2 
7 fc 



(3.4) 



are the PBCS number projection integrals. The 
PBCS correction term Ae k can be found in Ref. [29| . 
G(kkk'k'; 0), and F{kkk'k'\ 0) stand for the usual proton 
or neutron particle-particle (pp), and particle- hole (ph) 
matrix elements of the residual interaction V, i.e., 

G(klmn;J) = (kl; J\V\mn; J). 

F(klmn;J) = (fcr 1 ; Jl^lmn" 1 ; J). (3.5) 

Note that these relations are valid for both identical and 
non-identical particles. 

The forward-going (A M ), and backward-going (Y^) 
PQRPA amplitudes are obtained by solving the RPA 
equations 



A» B 
-B* -A*_ 



(3.6) 



with the PQRPA matrices defined as: 
A^pn^p'n'; J) 

= «JlWn' + N^{pn)N-^{p'n') 

x {[u p v n u p ,v n ,I z - 1+ »(pp')I N - 3 -»(nn') 

+v p u n v p ,u n ,I z - 3+ > 1 (pp r )I N - 1 -»{nn')}F(pn,p'n';J) 

+ [u p u n u p ,u n ,I z - 1+ »(pp / )I N - 1 -»(nn') 

+ v p v n v p ,v n ,I z - :i+ ^{pp l )I N - :i -^{nn')]G{p7i 1 p'n'- J)} 

B(pn,p'n'; J) 

1 I 2 (J„'\ rZ-lr l\ tN-2i 



= N-^{pn)NZ^{p'n<)r Z - 2 {pv')I N ~ 2 (nn>) 



x [(v p u n u p ,v n > + u p v n v p >u n >)F(pn,p'n; J) 
+(u p u n Vp>v n > + VpV n Up>u n >)G(pn,p'n'} J)]. 



where 



o = z-i+» 



• £ 



N-l-fi 



are the unperturbed energies, 

N^(pn) =I z - 1 +»(p)I N - 1 -i*(n) 
are the norms, 
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Rg(k) + RK(kk) 



I K (k) 



I K 



(3.7) 
(3.8) 
(3.9) 

(3.10) 



are the projected quasiparticle energies, and the quanti- 
ties R K are defined as [29J 



<(fc) 



Rf x {kk) 



^(2 Jfcl +l)< efcl ^- 2 (fcfc 1 ) 

ki 

1^(2^+1)^(2^+1)1/2 

kik2 

[vl l vl 2 F{k 1 k 1 k 2 k 2 ; 0)/^ 4 (A: 1 fc 2 fc) 
Uk 1 Vk 1 u k2 v k2 G(kik 1 k 2 k 2 ; 0)I K ^ 2 (k 1 k 2 k)] , 
e k [utl K (kk)-vll K - 2 (kk)} 

+ E ^+l^ K^ifci^O) 

[ul^^ihkk) - v 2 k I K - A {k 1 kk)] 

- Uk 1 Vk 1 UkV k G(kikikk;0)I K ^ 2 (kikk)} . 

(3.11) 

Both positive and the negative solutions are physically 
meaningful. For fi = ±1 the positive solutions describe 
excitations in the (Z±l, iV+ 1) nuclei, while the negative 
energy solutions represent the positive energy excitations 
in the (Z+ 1, A^± 1). Thus only one RPA equation has to 
be solved, either for [i = +1, or for /i = —1, to describe 
the excitations to the Z ± 1, iV + 1 nuclei. This is well 
known feature of the charge- exchange modes [47l - l50j . 

Let us be more specific, and take advantage of the in- 
dex / to label different final states | JJ) with same spin 
and parity. Evidently, / will run from 1 up the total 
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number / max of two-quasiparticle configurations \pnJ 7T ). 
Moreover, the eigenvalue problem (|3.6p has 2/ max solu- 
tions, and we will use the index F to label them. Thus, 
if fj, = +1 one have: 

• for u)+i(J F ) > (1 < F < / max ): 

X + i(pnJ f ) = X + i(pnJ F ), 
Y +1 (pnJ f ) = Y +1 (pnJ F ); (3.12) 

• for w+i(Jf) < (/ max < F < 2/ max ): 

= -w+i(Jf), 
X-i(pnJf) = Y+^pnJp), 
Y^{pnJf) = X* +1 (pnJ F ). (3.13) 

Finally, to store the eigenvalues and eigenfunctions it is 
convenient to define the index F as 



77 _ ) /' ^ — /max! /o i /) \ 

< 2/ max - / + 1, for F > / max . ^- i4 ^ 



B. QRPA 

The usual gap equations are obtained from Eqs. (I3.6[l - 
32) by: 

1. Making the replacement ejt — > — Ak, with Ak 
being the chemical potential, and taking the limit 
I K —> 1. That is, the Eq. (13.11) remains as it is, 
but instead of (13.21) and (13.31) one has now 



(3.15) 



and 



e k = ek - X k + ^2^±^-vl,F(kkk'k';0). 

(3.16) 

2. Impose the subsidiary conditions 

Z = J^(2j p + lfv 2 3p , N = ]T(2j„ + l) 2 vl , (3.17) 

ip in 

as the number of particles is not any more a good 
quantum number. 

In this way the usual BCS gap equations read 

2(e fc - X t )u k Vk = (wfc - v 2 k )A k . (3.18) 

This equation, together with the normalization condition 
Uf. +1*1 = 1, has as solution the occupation probabilities 
(for example, from the Chapter I of Rowe |47j ) 



i i 



efc - X k 
E k 



2 1 /. e fc - A fe 



which depend on the quasiparticle energies 

£ fc = ^ \e k - A fe ) 2 + A 2 k , (3.20) 
and the pairing gaps 

=-\Y, ^j^ wG(kkk'k'-,0). (3.21) 

The QRPA equations are recovered from (|3.6p by i) 
dropping the index n, ii) taking the limit I — > 1, and hi) 
substituting the unperturbed PBCS energies by the BCS 
energies relative to the Fermi level, defined by equation 



E^ — ±Ek + Ak, 



(3.22) 



where Ek are the usual BCS quasiparticle energies de- 
fined in p. 201) . In this way the unperturbed energies in 
(|3.7p are replaced by 



&=E jf +E Jn +n(\ p -X n ). 



(3.23) 



These energies, however, are not used in the QRPA eigen- 
value problem. Namely, the coefficients X(pnJf) and 
Y(pnJf), and the eigenvalues ui(Jf) are obtained from 



A B 
B A 



X 
-Y 



(3.24) 



where 

A(pnp'n';J) = (E p + E n )5 pp >S nn , 

+ (u p v n u p >v n > + v p u n v p 'U nl )F(pnp'n'; J) 
+ (u p u n u p >u n > + v p v n v p iv n >)G(pnp'n'; J), 

B(pnp'n';J) = (v p u n u p >v n > + u p v n v p >u nl )F(pnp'n'; J) 
+ (u p u n v p 'V n ' + v p v n u p > u n ')G(pnp' n '; J) . 

(3.25) 

In the pn-QRPA the eigenvalues occur in pairs ±o»(Jy), 
but the negative energies don't have a direct physical 
meaning. The perturbed energies for daughter [Z + 
fx, N — (i) nuclei are defined as 



J p (J f ) = u(Jf) + fi(X p - X n ). 



(3.26) 



There is, however, only one set of eigenfunctions 
(X(pnJf),Y(pnJf)) for both /! = 1, and fj, = —1. This 
is a very important difference in relation to the PQRPA 
case, which is crucial for the distribution of the transition 
strengths. 



C. Nuclear Matrix Elements 

When the excited states | J/) in the final (Z± 1, iV=p 1) 
nuclei are described within the PQRPA, the transi- 
tion amplitudes for the multipole charge-exchange op- 
erators (j2~2Tl) and (j2~2"2"j) read 



(J f ,Z- 



■/i,iV-/x||Oj||0H 



(/Z-l+ M (p)/«-l + M („))l/2 



(3.19) 



A_ M (pnJ) 



^Y*(pnJ f ) 



(3.27) 
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with the one-body matrix elements given by 

*„(!-- t )- <p j |0j||n) f ^ = +; , 

M ' v / 2TTT \ UnVp, for/i = -l 

(3.28) 

and J = J. 

In the QRPA case, using the limit I K -> 1 in (|3~2"71) . 
the nuclear matrix elements for the multipole charge- 
exchange operators Oj are 

(J/,Z + M,iV-/z||Oj||0+) 



(3.29) 



with the same one-body matrix elements (|3.28jl . 

The unperturbed and perturbed transition strengths 
are defined, respectively, as 



and 



Sl(pnJ) = \A^pnJ)\ 2 (3.30) 



S^Jf) = |(J/,Z + M,iV-Ml|Oj||0 + )| 2 (3.31) 



One might be particularly interested in the Gamow- Teller 
(GT) and Fermi (F) /3-decay strengths (B- values), in 
which case (|3 .30[) and p.31[) are evaluated for the op- 
erators Oj, which don't contain the radial form factors 
jj(p). That is, O = 1, and Oi = er, for F and GT oper- 
ators, respectively. We denote the B-values as S^pnJ) 
and S^(Jf). Occasionally one also might want to calcu- 
late the the energy distribution of the last one, i.e., 



S,(J f ,E) = IY, 



S^Jf) 



(3.32) 



where one usually takes 77 = 1 MeV. 



(b) conf.inc, which has the dimensions for the quasi- 
particle state configurations, the hamiltonian matrices 
(A, B) or (A, B) which are diagonalized, the forward and 
backward amplitudes, and the eigenvalues. 

There are two input files: 1) qrapout.dat, 
where is listed the output file that shows the neu- 
trino/antineutrino [vjv) cross section, as a function of 
the incident neutrino, or the muon capture rate, for each 
state of a given nuclear spin in the daughter nucleus, and 
2) the above mentioned qrapin.dat, which contains: a) 
the quantum numbers of all single-particle state (sps), 
and the corresponding single particle energies (s.p.e.), b) 
the mass and the proton number of the parent nucleus, 
c) the neutron and proton pairing strengths for the BCS 
approximation, d) the particle-particle, and particle-hole 
strengths of the residual interaction, e) the position of 
the Fermi level, and the experimental gap for neutrons 
and protons, and f) the Q-value for the vjv scattering. 

There are three default output files. Two of them, 
AUXI.OUT and OUT.OUT, contain the results of the 
nuclear structure model, whereas the results for the weak 
processes appear in the file created by qrapout.dat. 
For example, if one is interesting in the multipole 
Jj = 1 + with a single-particle space of six levels 
in 12 C ("set 1"), we can introduce in qrapout.dat 
the file names QNC.out (PNC. out) for neutrino 
capture, QAC.out (PAC.out) for antineutrino cap- 
ture, QMC.out (PMC. out) for muon capture, us- 
ing the QRPA (PQRPA) model. The auxiliary out- 
put files AUXI.OUT and OUT.OUT are relabeled to 
(QAUXI.OUT, QOUT.OUT) and (PAUXI.OUT, 
POUT.OUT) for QRPA and PQRPA respectively. 

All just mentioned outputs are included as examples. 
The following units are employed: i) 10~ 42 cm 2 , for neu- 
trino or antineutrino-nucleus cross sections, ii) 10 4 s _1 , 
for muon capture rates, and iii) MeV, for energies. 



A. Reading the data 



IV. COMPUTER PROGRAM AND USER'S 
MANUAL 

The QRAP code evaluates the electron neutrino- 
nucleus interaction described by equation (|2.1[) 
(IREAC=1 for NS, IREAC=2 for AS) and (|2~2"]) 
(IREAC=0 for MC). The processes, from the ground 
state of the even-even father nucleus (Z, N) to the 
excited states with spin (ISPIN) and parity (IPARI) 
in the odd-odd daughter nucleus [Z ± 1,N =p 1), are 
calculated by using the QRPA model (IQP=0) or the 
PQRPA (IQP=1). These options must be setup in the 
input data file, qrapin.dat, which is supplemented with 
two included files: 

(a) sp.inc, containing the dimensions of single-particle 
quantum numbers, occupation probabilities, quasiparti- 
cle quantities and strength amplitudes for allowed tran- 
sitions; 



There are three sets of input data in qrapin.dat sep- 
arated in modules labeled as: *Data set 1 for a single- 
particle space of six levels in 12 C (0, 1, and 2 ftjjj oscillator 
shells) , *Data set 2 for a single-particle space of ten lev- 
els in 12 C (0,1,2, and 3 Huj oscillator shells), and *Data 
set 3 for single-particle space of 12 levels in 56 Fe (2,3, 
and 4 hu oscillator shells). 

For each one of these input data, the number of sps rep- 
resents the available space where one wants to solve the 
BCS (or PBCS) problem given by equations (13. 17)) and 
(I3.18P f (|3.ip - (l3.2p ). It contains the necessary number of 
harmonic oscillator shells leading to a smooth smearing of 
the Fermi's surface. The Fermi level with the neighboring 
levels constitute the active shell for the mentioned smear- 
ing. For example, in 12 C (ground state with J = + ) the 
active shell is composed by the lp^/2 and lpi/2 levels. 
According to the single-particle shell model the sps filled 
up the IP3/2 orbital, and the nucleons can be promoted 
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to the lpi/2) creating a particle state in lpi/2 and a hole 
state in lp^/2- This scheme describes the first particle- 
hole (ph) excitation on 12 C in order to obtain the 12 N or 
12 B ground state with J = 1 + , by promoting a proton or 
a neutron, respectively. 

Let us show, as example, a data input of six sps for 
12 C: *Data set 1. The rows starting with a symbol "*" 
are not read as input and just serve to remind the user 
on the meaning of the physical quantities. Taking out 
the comments "*" in the first lines of this file, we have 

1 +1 0.0 1 1 
06 06 1 1 
101 -20.09 112 -6.02 111 -0.29 123 3.07 201 3.85 122 7.18 
101 -18.19 112 -3.17 111 2.79 123 5.73 201 6.06 122 9.36 
12 06 28.80 28.85 30.0 50.0 27.0 64.0 
2 1 6 6 1.00 6.88 

2 1 6 6 1.00 7.00 

17.338 

First line: Nuclear spin (ISPIN=1) of the daugh- 
ter nucleus, parity (IPARI=+1), coupling strength of 
particle-particle channel (t = 0.0) (see definition below), 
index of neutrino reaction (IREAC=1) and the index 
(IQP=1) to solve the PQRPA problem. 

Second line: Number of neutron sps (NSN=06), 
number of proton sps (NSP=06), index to solve the 
QRPA equation in the Tamm-Damcoff approximation 
(ITD: no, 1 yes), index to print the matrix elements 
of the nuclear hamiltonian to be diagonalized (MAPR: 
no, 1 yes), index to solve the BCS equation with the 
self-energy term (IFMU: no, 1 yes), index to make the 
QRPA matrix with the branch (Z + 1, N - 1) (fi = 1) or 
(Z - l,N + 1) Git = -1) of PQRPA (IPRO: (/x = 1), 1 

(m = 7 i)). 

Third and Fourth lines: quantum numbers and the 
s.p.e. for each neutron and proton sps, respectively. 
They are represented in the same way as in the shell 
model scheme, with their respective quantum numbers 
(n I (j + 5))- For example, 101 — > lsi/2 where 1 is princi- 
pal quantum number, corresponding to the first harmonic 
oscillator level (n), corresponds to the orbital angular 
momentum I = s and the last number l=j + i = i + i. 
Table U shows the notation and the corresponding quan- 
tum numbers, as well as the PBCS quasiparticle energies 
ef , and ef , defined by (|3TTUj) . 

Fifth line: Mass number A (IAM=12), proton 
number Z (IZ=6), and the following coupling con- 
stants: 1) neutrons and proton pairing v pairN (Vs- 
pairN=28.80), and v pairP (VspairP=28.85), 2) singlet 
and triplet particle-particle (pp) v pp (VsPP=30), and v PP 
(VtPP=50), and iii) singlet and triplet particle- hole (ph) 
v ph (VsPH=27), v ph (VtPH=64). 

Sixth and Seventh lines: Position of the Fermi level 
(LEVEL=2), initial (IIQ=1) and final (IFQ=6) states 
for which the BCS equations must be solved, number 
of particles interacting (NPIQ=6) in these levels, and 
the experimental gap (DELTAQ=6.88 or 7.00) defined 
below in equation (|4.2p ; for neutrons and protons, fifth 



TABLE I: Notation for the quantum numbers, the resulting 
quasipartice energies ef for neutrons and ef for protons, and 
the pairing strength within the PBCS. The energies are 
given in units of MeV, and v® mr is dimensionless. 



Notation 


Shell 


n 


£ 


j + 1/2 


c j 


ef 


101 


1*1/2 


1 





1 


-20.09 


-18.19 


112 


1^3/2 


1 


1 


2 


-6.02 


-3.17 


111 


lPl/2 


1 


1 


1 


-0.29 


2.79 


123 


W5/2 


1 


2 


3 


3.07 


5.73 


201 


2S1/2 


1 





1 


3.85 


6.06 


122 


W3/2 


1 


2 


2 


7.18 


9.36 


pair 


28.80 


28.85 



TABLE II: Spin and parity for the one-quasiparticle space 
used in the input for 12 C. 





lSl/2 


1P3/2 


lfl/2 


2S1/2 


lds/2 1^3/2 


lSl/2 


0+1+ 










1P3/2 


l-,2" 


0+1+ 
2+,3+ 








lPl/2 


o-,i- 


1+2+ 


0+1+ 






2S1/2 


0+,l+ 


l-,2" 


o-,i- 


0+1+ 




lds/2 


2+,3+ 


l-,2" 


2",3" 


2+,3+ 


0+,l+,2+ 




3", 4- 






3+,4+,5+ 


1^3/2 


l+,2+ 


o-,i- 


l-,2" 


l+,2+ 


l+,2+ 0+,l+ 




2-,3" 






3+,4+ 2+,3+ 



and sixth lines respectively. 

Eighth line: Q- value minus the lepton mass for v jv 
scattering (EGS=17.338 for 12 N jp). It can be fixed as 
being the energy of the ground state in the daughter nu- 
cleus. The lepton mass must be added to EGS to obtain 
the Q-value for the reaction. 

B. Running the code 

As first step the QRAP solves the BCS problem. In 
this case, one needs to adjust the pairing strength to 
reproduce the experimental pairing gap. 

Next, one can solve the PBCS problem or directly 
the QRPA if the option IQP=0 was selected. If IQP=1 
then the PQRPA equations are solved. It means that 
QRAP firstly calculates the nuclear matrix elements in 
the QRPA or PQRPA by selecting the option IQP=0 or 
1, appropriately. The option for which type of weak in- 
teraction process one wants to evaluate is adopted with 
IREAC in the input data. We recommend first to adjust 
the pairing strength as it is explained below. After this it 
is convenient to fit the parameters of the residual inter- 
action using the option IREAC=3 for the muon capture 
rate because this calculation is fast. Physically you can 
check quickly how good is your choice of parameters be- 



cause the values for inclusive muon-capture rate, and GT 
B-values are available in the literature (see for example 
Refs. HHH). 

For the residual interaction the code assumes a delta 
force, 



V 



-4tt (v a P s + VtP t ) S(r), 



(4.1) 



which has been used extensively in the literature [26 
to describe single and double beta decays. 

Next, we explain how the parameters of the interaction 
are adjusted using for example the input data for six 
levels in 12 C. The results are presented in output file 
OUT.OUT. 

• Adjusting the gap A&: 

The parameters yP a ' lrN and v pa%rF are adjusted to re- 
produce the experimental gap A N for neutrons, and A z 
for protons, by solving the BCS equations (|3.17p and 
(|3.18[) in a self-consistent way. The experimental gaps, 
according [13, Eq. 2.96], are: 



A N = --{B(Z, N —1) — 2B(Z, N) + B(Z, N + 1)} . 



{B(Z - 1, N) - 2B(Z, N) + B(Z + 1, N)} , 



(4.2) 



where B(N, Z) is the binding energy of the even-even nu- 
cleus (Z, N). This is the most common fit (Fitl) used in 
several works with standard QRPA 28-30]. In this case, 
the A N ( Z ) must be equal or approximately equal to the 
energy A^^ FL of the corresponding to the Fermi level 
(FL). To solve the set of PBCS coupled equations (l3~Tj) - 
Q3.4[) for Uk and Vk it is recommended to obtain first the 
solutions for the BCS problem, as these probability occu- 
pations are use as input for the PBCS case. The PBCS 
coupled nonlinear equations are solved consistentl y w ith 
Powell Hybrid method using subroutine HYBRD [55]. 

The results of the BCS or PBCS problem are shown as 
tables in the first lines of OUT.OUT for neutrons and 
protons, respectively. The quantities defined by (|3.10p 
and (|3.11[) are presented there. In particular, the pro- 
jected quasiparticle energy defined in (13. 10)) are 



PROYSP 



E(+) 
E(-) 



_K-2 



with k above Fermi level, 
with k below Fermi level, 

(4.3) 



which means that E(+) corresponds to a particle state, 
and E(— ) to a hole state. The values of A^ are shown 
in the ninth column of the table labeled as CONFIGU- 
RATION SPACE. This Fitl comes from the fact that the 
experimental energy difference between the states that lie 
just above (p state) and just above (h state) the FL is 
approximately twice the experimental gap, i.e., 



E 



K 



2A 



K 



(4.4) 




FIG. 1: Ground state energy E, and B(GT)-value in 12 N for 
different couplings r in the p/i-channel, as a function of the 
pp-channel coupling t. The experimental data [5ll. |58|. [5^] are 
also shown. The value r — 2 corresponds to v^ h = 27 a nd 



= 64 (Case PII of [23||). The Data set 2 was employed. 



There is another fitting procedure for the pairing gap 



that is called by Fit2. In Fit2, all the s.p.e. e 



N(Z) 



from 



Table I are varied with a % 2 search to account for the 
experimental spectra Ef 



Z(N) 



E z{ - N) = E(+), for a particle state, 



Z-2(N-2) 



E 



Z(N) _ 



= E(— ), for a hole state. 



for K = N or Z. 



In Fit 2, the Eq. (|4.4[) is automatically satisfied. This 
procedure was employed to obtain the ej spectra shown 
in [H, Table III], whereas the ej for the reduced space 
of six levels in the present example are shown in Table HI 
These s.p.e. are used in input data Data set 1. 

To make the calculations as simple as possible the Fitl 
procedure is the usual choice, with the ej spectra ob- 
tained either from a harmonic oscillator or from a Wood- 
Saxon potential, and by varying the coupling vP airN and 
v pairP tQ satigfy the condition A^ ] FL w A N W. 

For 56 Fe the input data is called *Data set 3. The s.p.e 
of the active 3ftu; shell were taken from the experimen- 
tal energies of 56 Ni, and the 2hw and Afku shell energies 
were taken from the harmonic oscillator energies with 
hu/MeV = 45 A 1 / 3 - 25 ,4 2/3 . Fitl was employed to 
adjust the experimental A N ^ for 56 Fe. 

• Adjusting the particle-hole couplings r and p 
In the particle-hole matrix element F, defined in 
Eq. (|3.5I) . the couplings v s and vt appear as linear com- 
binations v s + Vt and 3v t — v s . Therefore, it is convenient 
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FIG. 2: Results from PQRPA calculations obtained in 
Ref. [13], as a function of the pp parameter t, compared with 
the experimental data taken from Refs. [5lL l58l. [59|] . for: (i) 
ground state energy in 12 N (upper panel), (ii) _B(GT)-value 
for the f3 transition in 12 N (middle panel), and (iii) exclusive 
muon capture rate A(l]*") in 12 B (lower panel). The parame- 
ters of the ph channel for the 5-interaction are = 27 and 
64. 



ph 



to introduce the dimensionless parameters 



For nuclei with N > Z the pp-couplings are fixed on 
the basis of the SU(4) and isospin symmetry, as v PP = 
v pair , and v PP > vPP (^,[19. However, in Ref (H it was 
shown that this parametrization might not be suitable for 
N = Z. In fact, the best agreement with data in 12 C was 
obtained when the pp-channel is totally switched off, i.e., 
v PP = v PP = 0, and three different set of values for the 
ph-coupling strengths were used. These conditions are 
related with s = t for 12 C (N = Z), and with s = 1 and t 
variable in nuclei with N > Z. For 56 Fe were adopted the 
values s = 1 and t = 0. In the code QRAP the following 
conditions are standard: (i) s = t with t as a variable 
parameter for N = Z; and (ii) s = 1, and t as a variable 
parameter for N > Z, i.e., the residual interaction is 
defined as a function of two adjustable parameters v ph 
and t. 

Several experimental data are available in the liter- 
ature that can be used for fixing the residual interac- 
tion coupling constants, such as: ground state energies 
of daughter nuclei, £?(GT)-values for the 5+ or (3~ decay, 
and partial muon capture rate [SJ [H, [59| . 

One can use the reduced space of six levels to identify 
in the output file, the quantities shown in Figure [1] The 
results for three values of t are shown in Table IIIII 

TABLE III: Evolution of ground state energy, B{GT) and 
exclusive muon capture rate in 12 C, as function of the pp- 
channel parameter t. With [a] and [b] we denote, respectively, 
the output files "AUXI.OUT" and "PMC.out", with co fl (l+) 
is in units of MeV, Sp(Vf) is dimensionless and (A(lf), A) is 
in units of 10 4 s _1 . The parameters for the p/i-channel are: 



,,ph 



+ v. 



ph 



3vf h 



„ph 



2v, 



pair 



2 v 



pair 



where 



V 



pairN 



■ V 



pairP 



(4.5) 
(4.6) 

Moreover, to start we can use v ph = 27 and v ph = 64 (in 
units of MeV.fm 3 ). These values are inferred from the 
systematic study of energetics of the GT resonances done 
by Nakayama et al. [5(| (see also Ref. [3), and have been 
extensively used in the QRPA calculations of the /?/?- 
decay in 48 Ca [1^, [3(J H3|- Moreover, it makes sense to 
take the singlet ph coupling to be equal to v pair , obtained 
from the proton and neutron gap equations, i.e., 



y ph = v pair_ 



(4.7) 



• Adjusting the particle-particle couplings s and t 
Here also it is convenient to normalize to v pair the cou- 
pling constants v pp and v PP that appear in the pp matrix 



elements G in Eq. (|3.5p . and correspond, respectively to 
the channels (T = 1, S = 0) and (S = 1,T = 0), i.e., 



,J'i' 



pair 



t = 



pair 
Vs 



(4.8) 





State [File] 


Observable 


t 

0.0 0.3 0.6 Exp. 


12 N 


16 [a] 


w +l(!l6) 


18.319 17.951 14.970 17.34 [51] 




16 [a] 


S+i(lt 6 ) 


0.496 0.696 0.840 0.466 [58] 




1 [a] 


w-i(lfe) 


12.528 12.126 9.202 13.36 [51] 




1 [a] 




0.502 0.693 0.837 0.526 [58] 




16 [b] 


A+iUie) 


0.689 0.936 1.119 0.62(3) [59] 




[b] 


E f A+i(l+) 


1.722 1.537 1.183 



The values of and 5^(lt) in 12 N and 12 B can 

be found in the output file AUXI.OUT. In the present 
case the largest value of index / is / TOO a:=16. Both set of 
states, with /i = +1, and fi = — 1, are ordered from high- 
est to lowest energies. In the PQRPA, the most collective 
ones are that of the corresponding ground states: 1 1 J_ 16 ) 
in 12 N (and |l£_ 17 ) in 12 B) although there also are signif- 
icant strengths in the states F = 7, 11, and 14. In QRPA, 
the ground state is in 1 1 £ =16 ) for both 12 N and 12 B. These 
wave functions are presented below. For the PQRPA 
case, we also show the unperturbed energies w°(pnl + ) 
(which are not ordered), and the corresponding single- 
particle GT strengths S®(pnl + ), given respectively, by 
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flOJ) , and (f3T30|) for the GT operator er. The largest 
ones are S° +1 (lpl /2 , lp» /2 , 1+), and S!^(lpJ /2) lp» /2 , 1+), 
which in the particle-hole limit correspond to excitations 
1^3 , 2 — > lPi/21 an d IP3/2 ~~ > lPi/2- For s P ms an d parities 
7^ 1 + , or + , in AUXI.OUT are shown the energies 
ojft, but not the strengths S^. 

The results for the eigenvalue problem are displayed in 
the output file OUT.OUT. For the option M APR = 1 
are printed out the matrix elements (A, B) Eq. (|3.25[) for 
the QRPA, or (A^ =1 ,B) eq. fl3J]) for PQRPA. The nu- 
clear wave functions (X(pn; Jf), Y(pn; Jf)) are grouped 
to four, with the index F, defined in (|3.14[) , going from 
1 to f m ax in the QRPA case, and from 1 to 2f rnax in the 
PQRPA case. To make easy reading together with each 
set of wave functions are also printed: the value of /, 
the two quasiparticle configurations {p and n), and the 
unperturbed and perturbed energies. 

Recalling Eqs. flU), (|3~T2]) . and (j3~13l) for the energies, 
one discovers without difficulty that within PQRPA: 

1) The ground state in 12 N, with energy w+i(l + ) = 
18.319 MeV, has / — F = 16, and that its wave function 
is: 

| 12 N) = 0.963|l^ /2 l^ /2 ,l+)+0.232|l^ /2 l^ /2 ,l+) 
+ 0.122|1^ /2 1^ /2 ,1+) +0.105|lp 1/2 lK/ 2 ,l + ) 
+ ••• (4.9) 

2) The ground state in 12 B, with energy w_i(l+) = 
12.528 MeV, has / = 16, F = 17, and that its wave 
function is: 

| 12 B) = -0.971|lp^ /2 lp^ /2 ,l+>+0.204|l^ /2 l^ /2 ,l+) 
- 0.125|1^ /2 1K /2 ,1+) + 0.090|K /2 lp 1/2 ,l+> 
+ ••• (4.10) 

One proceeds in a similar way for the QRPA output, 
with energies now given by Eqs. (|3.23[) and (|3.24[) . Now, 
the ground state energies in 12 B, and 12 N, are, respec- 
tively, = 12.437 MeV, and w+i(l+) = 17.992 
MeV, while the wave function for both nuclei is: 

|lf 6 > = -0.272|lpJ /2 lpy /a ,l+)-0.759|l^ /2 l^ /2) l+> 
+ 0.356| lp v 1/2 lp v 3/2 , 1+) - 0.472|l^ /2 l^ /2 , 1+) 
+ ••• (4.11) 

From the comparison of the wave functions (|4.9p . and 
(|4.11[) it can be easily figure out why Volpe et al. [2lJ 
called attention to "difficulties in choosing the ground 
state of 12 N, because the lowest state is not the most 
collective one" when the QRPA is used. This is an 
important issue that clearly gives you an idea about 
the need for the number projection. In fact, as seen 
from fll?)) . and P~TTj) . the PQRPA yields the correct 
one-particle-one-hole (lplh) limits lpy 2 — > IPin an d 

1^3/2 ~~ * ^1/2' ^ or 12 ^' ano ~ ground states, re- 
spectively. All remaining configurations comes from the 



higher order 2p2h, and 3p3h excitations. Contrary, the 
QRPA state (|4.1ip is dominantly the two-hole excitation 
[(±Pg/ 2 ) _1 , (lp^ 2 ) -1 ], which corresponds to the ground 

state of 10 B. This should not be a surprise, as we know 
that the proton-neutron QRPA states are the same for 
all four nuclei 12 N, 10 B, 14 N, and 12 B. More details on 
this question can be found in [23], Fig. 3]. The lplh 
amplitudes[(l^ /2 )" 1 , lp" /2 ], and [(lp£ /2 )~\ (lPi/ 2 )] are 
dominantly present in the QRPA states / = 13, and 
/ = 15, i.e., 

|1+) - -0.476|l^ /2 K /2 ,l+)+ 0.437|lpJ /2 l^ /2 ,l+) 

+ 0.441|lpJ /2 lp£ /2 , 1+) - 0.096|lpJ /a l# /a , 1+) 
+ ••• 

|1+) = 0.703|lpJ /2 l P r /2 ,l + )+0.708|l^ /2 l^ /2 ,l+) 
+ ••• (4.12) 

The wave functions displayed above clearly evidence 
the superiority of the PQRPA on the QRPA. 

• Output for the j/- nucleus processes 

The output of the results for the weak processes is 
selected according to the value of IRE AC: 

IREAC=0 prints the results for the muon capture 
rate in the file (QMC.out or PMC.out). For J* = 0+ 
or J n = 1 + are shown in in this output file the folded 
strengths S^(Jf,E) (S" TILDE) defined by ((3~52l . where 
'ENERGY' represents E. The partial capture rate for 
each state /, the perturbed energy w^ = _i, and the 
strength S M =_i(J/) (if J w = 0+, or J T = 1+) are shown 
in the table labeled CAPTURE RATE. The total capture 
for the evaluated spin J*' is presented in the last line. 

IREAC=1 or IREAC=2 prints the results for 
the neutrino or antineutrino cross sections in the files 
(QNC.out/PNC.out) or (QAC.out/PAC.out). We repeat 
in this output the folded strengths 5 M (J/,S) (|3.32j) for 
J n = 0+ or J n = 1+. The cross sections (SIGMA(Enu)) 
are calculated as a function of the neutrino energy (Enu) 
for each nuclear spin from / = 1 to / = / max . The per- 
turbed u}/j,=±i energies for the daughter nucleus are also 
shown according the process related. The absolute value 
of maximum (cos 9 = — 1) and minimum (cos 9 = 1) nu- 
clear momentum transfer (|k|) in units of MeV/c for each 
energy are also printed. 

Note: The cross sections are printed up to a maximum 
energy of 250 MeV. Depending upon the single particle 
space employed, the cross sections, as a function of the 
neutrino energy, should be restricted to lower energies. 
This issue will be discussed and explained in details in 
a next work [6(|. Anyway, the PQRPA cross sections 
obtained within the single particle space provided as ex- 
amples are well behaved up to E v /a < 100 MeV on av- 
eraged according to *Data set 1, *Data set 2 and *Data 
set 3. This interval of energies is important for supernova 
neutrinos and low-energy decay-at-rest neutrinos [loj . 
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V. ROUTINES INCLUDED WITH THE CODE 

QRPA solves the pn-QRPA or pn-PQRPA charge- 
exchange problem for a nuclear spin JJ of the daughter 
odd-odd nucleus. 

SUAVE calculates and prints the folded strength 
,r = 0+ or r = 1+ given by Eq. (j3~32l) . folding the 
Sn(Jf) stenght with a Lorentzian function with rj = 1 
MeV. 

RMUONCAP calculates the muon capture rate 
given by formula (|2.29p . 

SIMPSN2 calculates the neutrino or antineutrino 
cross sections as a function of the neutrino energy. This 
subroutine uses the function G to call the subroutine 
SECCION, which evaluates the cross section formula 
(|2.23[) using the Gauss-Legendre N-point quadrature for- 
mula [61( on the function F to evaluate the angular inte- 
gration of the transition amplitude times Eg. 

MATRIXP computes the matrix elements with the 
delta residual interaction given in Ref. [62| for the 
PQRPA. The matrix elements were modified according 
to the projection procedure shown in Eq. (|3.7[) . 

MATRIX computes the matrix elements with the 
delta residual interaction given in Ref. [62| for the QRPA. 
The matrix elements are shown in Eqs. (|3.25l) . 

RPA finds eigenvalues and eigenvectors for the QRPA 
or PQRPA equations. It uses the subroutine EIGRF 
and other related subroutines from the IMSL Library [(33| 
to orthonormalize the eigenvectors. 

GAPII solves the set of BCS coupled Eqs. (|3TH) and 
(|3.18p to obtain the and Uk for neutrons and protons. 

CONFGT builds up the pn configurations for a given 
spin and parity. 

FAUX evaluates the particle-particle matrix elements 
G(kkk'k'; 0) and F(kkk'k'; 0), which are used to solve the 
gap equations for neutrons and protons using the delta 
interaction. 

RADWF computes harmonic oscillator radial wave 
functions. It uses the additional subroutine OSCILL to 
evaluate the radial coefficients. 

HYBRD finds a zero of a system of N nonlinear equa- 
tions in N variables by a modification of the Powell Hy- 
brid method. This subroutine was provided by the Ar- 
gonne National Laboratory 55]. It uses the subroutine 
FCN to calculate the PBCS nonlinear equations given 
by formula (13. ip . 

FKPERMAT evaluates the perturbed matrix el- 
ements for the weak decay operator, according to 
Eq. (ET2T|) for PQRPA, and Eq. (ET29)) for QRPA. The 
radial part of the SPNME were defined in Ref. [ill. This 
subroutine uses the subroutine ANGULARMATRIX 
to calculate the angular part of single-particle matrix el- 
ements defined in Ref. [57] . and shown in the Appendix 
A for the sake of completeness. 

There are other routines in the code that are shortly 
described as follows. PRINMA prints the matrix el- 
ements (A,B) for the QRPA, or (A^B) for PQRPA, 
SKIPCOM is used to skip comments in the input file, 



UNPMOM3 evaluates the unperturbed projected ma- 
trix elements for beta decay, BETMAT2 is used to cal- 
culate the single-particle matrix elements for beta decay, 
PROENER calculates the quantities for the projected 
quasiparticles energies in (|3.10[) . 



VI. THINGS TO DO 

1. Use the sample input Data set 1 to obtain the results 
presented in Table ITTTT 

2. Modify the input Data set 1 by Data set 2, set- 
ting all parameters of the residual interaction to zero. 
These values correspond to BCS or PBCS approxima- 
tion. Compare the folded strength of Data set 1 with 
Data set 2 shown in Fig. 4 of Ref. [HI. 

3. In Ref. [65[ the s.p.e. for neutrons were changed to 
analyze the systematics of the paring strength in the odd 
carbon isotopes. Change the s.p.e for neutrons in Data 
set 2 and reproduce the systematics shown in the level 
scheme of Fig. 2 and the spectroscopic factors of Fig. 3 
of Ref. [H. 

4. Compare the QRPA and PQRPA results for the ex- 
clusive v e — 12 C cross section, as a function of the neutrino 
energy, with the DAR experimental data from Ref. [68j . 
Note that the QRPA result is not collective, and the ad- 
dition of other 1 + cross sections (for example, that of 
states l|4.12p ) is required to get agreement with the ex- 
perimental value. 
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Appendices 
A. Single-particle nuclear matrix elements 

The elementary operators defined in Eq. (|2.16l) have 
the reduced single-particle pn matrix elements (RSPME) 
defined in [4l|, [HjJ (Recall that k = |k| and p = — iV). 

For the RSMPE dependent on the tensor product of 
spherical harmonic times the nucleon velocity we have 

(P,Qp |UplliL(Kr)[YL(r)<g> V]j||n,(J„ ±),j„) = 
( ~ 1 J^ F ' +L W[j } {pn)R ( ~ ] {pn; n) + {pn)R ( L +) (pn; k) , 

(6.1) 
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with angular and radial parts, respectively: 



2 ^ T 



( /„L|/„ \ '*■ \ 



R^ipn;*,) = / M„ p ,z p (r) 
Jo 



± 



L J 1 

In In ~F 1 

2/„ + 1 ± 1 



2r 



(6.2) 



We use here the angular coupling J = V2J + 1 

and (l p L\l n =F 1) is the short notation for the Clebsh- 
Gordon coefficient (l p 0L0\(l n =p 1) 0). 

For the scalar product of spin times nucleon velocity, 
we have 

(P,(1 P |)J P ||jj(/cr)Yj(r)(<r- V)||n,(J n ±),j„) =(6.3) 



i 

/47T 



(pn)i^ } (pn; k) + Wj (+) {pn)R { j +) {pn; k) 



?(-)/ 



with the angular part 
^ (±) (pn) = ±(-l) i "+^+ J + 1 /2V6jf p j pJ „(/ rl + | T i) 1/2 



1 



T 1 jn 2 I (6.4) 
jp ^p J 



being the radial part i?j (pn; ft) as in (|6.2p . 

The RSPME of the two operator independent of the 
nucleon velocity are written below. For the the spherical 
harmonic operator we have 

(P, (l P |) J P ||jj(/tr)y,(r)||n, (f„ i),j„> = (6.5) 
with the angular and radial parts, respectively: 



W M ( P n) = (-i)^- J "Jjpj„ [' J f ^ „ 

2 2 U 



(6.6) 



iq(pn; k) 



u n P ,ip( r ') u ™„j„(r) jj(Kr) r 2 dr. 



Finally, the RSMPE dependent of the tensor product of 
spherical harmonic times the spin operator reads 



(P> (lp §), JpIIjl(w) [F L (r) ® crjj 1 1 7i, (/„ |),j n ) 



(-i) L 



-W L j(pn)i?p(pn;fi;), 



V4ir 

where the angular part is 

W u (pn) = (-l) /f, V / 6jpj„fpf„j p L J 

2 ^p ip 

- Z 7 

1 L J 



(6.7) 







(6.8) 



with the radial part i?"(pn; k) given by (|6.6 



B. Fermi function and effective momentum 
approximation (EMA) 



To account for the Coulomb interaction between the 
charged lepton and the residual nucleus, the QRAP code 
is setup to use by default the Fermi function [45], 46] . This 
correction was employed in several works for reactions on 
12 C with neutrinos from the DAR of p + . As pointed out 
in Ref. [2l|, the quantity piRa is of the order of 0.5, 
where pi is the lepton momentum, and Ra is the radius 
of the nucleus. Thus, the correction is well described by 
a Fermi function. Yet, for high energy neutrinos, e.g. 
neutrinos from the DIF of ir + , the outgoing muons have 
Ps.Ra > 0.5. For these relativistic leptons, the effective 
momentum approximation (EMA) [661 ] should take care 
of the Coulomb field of the daughter nucleus, instead 
of the Fermi function. This prescription for the Coulomb 
correction is considered in the code within the subroutine 
SECCION. More precisely, with EMA=0 the Fermi func- 
tion is employed, while with EMA=1 the EMA prescrip- 
tion is used. In the EMA procedure, the lepton energy 
and momentum are modified by a constant electrostatic 
potential within the nucleus 



Ei,eff = E £ ~ Veff, 



Pt,eff = \l E f 



t,eff 



with Veff = 41/ c (0)/5 = -6Z f a/5R A [3J, ISJj- These 
two approximations for the Coulomb correction were 
tested in the calculation of the inclusive cross section for 



neutrino scattering on U8 Pb 



As shown in Ref. 1 38 



the Fermi function correction overestimates the cross sec- 
tions at higher neutrino energies where the EMA provides 
a more reliable approach. Thus, we recommend to use 
the Fermi function correction in the range of neutrino en- 
ergies for which the cross section is below the correspond- 
ing EMA value, whereas the EMA could be employed at 
higher energies, as shown in previous studies [2l[ |38|. 

As a final comment, the QRPA code could be easily ex- 
tended to calculate ^-induced processes. This was done 
in Refs. [22|, |23j to calculate v^ — 12 C cross sections using 
the EMA prescription for the DIF regime of the LSND 
experiment. The nuclear structure calculations remain 
the same, while the kinematics changes by changing the 
electron mass to the muon mass in the variable RMLEP. 
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